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In  the  context  of  stabilization  of  high  order  spectral  elements,  we  introduce  a  dissipative  scheme 
based  on  the  solution  of  the  compressible  Euler  equations  that  are  regularized  through  the  addi¬ 
tion  of  a  residual-based  stress  tensor.  Because  this  stress  tensor  is  proportional  to  the  residual 
of  the  unperturbed  equations,  its  effect  is  close  to  none  where  the  solution  is  sufficiently  smooth, 
whereas  it  increases  elsewhere.  This  paper  represents  a  first  extension  of  the  work  by  Nazarov 
and  Hoffman  [Nazarov  M.  and  Hoffman  J.  (2013),  Int.  J.  Numer.  Methods  Fluids,  71:339-357 .] 
to  high-order  spectral  elements  in  the  context  of  low  Mach  number  atmospheric  dynamics.  The 
simulations  show  that  the  method  is  reliable  and  robust  for  problems  with  important  stratification 
and  thermal  processes  such  as  the  case  of  moist  convection.  The  results  are  partially  compared 
against  a  Smagorinsky  solution.  With  this  work  we  mean  to  make  a  step  forward  in  the  implemen¬ 
tation  of  a  stabilized,  high  order,  spectral  element  large  eddy  simulation  (LES)  model  within  the 
Nonhydrostatic  Unified  Model  of  the  Atmosphere,  NUMA. 


1  Introduction 

Recently  [17],  a  numerically  stable  and  computationally  inexpensive  large-eddy  simulation  (LES) 
model  for  compressible  flows  was  designed  for  adaptive  finite  elements.  It  is  a  close  relative  of  the 
entropy- viscosity  method  by  Guermond  and  co-workers  (see,  e.g.  [7]),  although  no  entropy  equation 
is  used  to  construct  the  dynamic  viscosity  coefficient  of  the  stress  tensor. 

In  the  current  paper,  we  explore  the  capabilities  of  the  aforementioned  LES  model  to  act  as  a  sta¬ 
bilization  method  for  the  spectral  element  solution  of  the  Euler  equations  at  the  low  Mach  number 
regimes  typical  of  atmospheric  flows.  This  effort  is  justified  by  the  fact  that,  within  the  community 
of  atmospheric  modelers,  there  is  still  a  widespread  concern  about  the  most  proper  stabilization 
scheme  to  be  used  with  either  Galerkin  or  other  approximation  methods  of  the  equations  of  atmo¬ 
spheric  dynamics.  Although  the  use  of  residual-based  stabilizing  schemes  has  been  largely  assessed 
for  the  finite  element  method  during  the  past  thirty  years  (e.g.  Streamline-Upwind/Petrov-Galerkin 
(SUPG)  [3],  Galerkin/Least-Squares  (GLS)  [10],  Variational  Multiscale  (VMS)  [9,  8,  2,  14]),  hy¬ 
per  viscosity  is  still  today  the  most  classical  approach  in  spite  of  its  important  drawbacks  and 
mathematical  inconsistency. 

This  work  is  a  first  step  toward  the  implementation  of  a  stabilized  high  order  spectral  element 
LES  model  (LES-SEM)  for  the  Nonhydro  static  Unified  Model  of  the  Atmosphere  (NUMA)  developed 
by  the  authors  [11,  5].  The  rest  of  the  paper  is  organized  as  follows.  The  set  of  equations  and  the 
LES  model  are  described  in  Section  2.  Some  basics  on  the  space  and  time  discretization  of  these 
equations  is  reported  in  Section  3,  which  is  followed  by  the  numerical  tests  and  results  in  Section 
4.  Some  conclusions  are  given  in  Section  5. 
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2  Equations  for  wet  dynamics 

Let  D  E  M3  be  a  fixed  three  dimensional  domain  with  boundary  dQ  and  Cartesian  coordinates 
x  =  (, x ,  y,  z).  Let  us  identify  the  dry  air  density,  the  velocity  vector,  and  the  potential  temperature 
with  the  symbols  p,  u,  and  8.  Let  us  also  define  the  mixing  ratios  of  water  vapor,  cloud  water, 
and  rain  as  qv  =  pv/p,Qc  =  Pci P  and  qr  =  pr/p,  where  pv,c,r  are  the  densities  of  vapor,  cloud, 
and  rain.  Furthermore,  let  p'(i,x)  =  p(f,  x)  —  po(z),  d'(t,x )  =  0(t,x)  —  6o(z),  and  p'(t,  x)  = 
p(t,x)  —  po(z)  be  the  perturbations  of  density,  potential  temperature,  and  pressure  with  respect  to 
a  hydrostatically  balanced  background  state  indicated  by  the  subscript  0.  Then,  the  strong  form 
of  the  time-dependent  Euler  equations  with  gravity,  g,  can  be  written  as: 


p't  +  u  •  Vp  +  />V  •  u 
u*  +  u  •  Vu  +  ±  V  •  (I pt) 
O'f  +  u-  V6» 


0, 

g{ i  +  tq-u  -  qc  -  Qr) k, 
Se(p,6,qv,qc,qr), 


(1) 


qit  +  u-X7qi  =  Sqi(p,8,qv,qc,qr),  for  i  =  v,c,r, 

where  I  is  the  identity  matrix,  k  is  the  unit  vector  [0  01]7”,  and  e  =  R/Rv  is  the  ratio  of  the 
gas  constant  of  dry  air,  R  and  the  constant  of  water  vapor,  Rv.  Because  moist  air  contributes 
to  the  buoyancy  of  the  flow,  the  right  hand  side  of  the  momentum  equation  is  corrected  with 
total  buoyancy  B  =  <7(1  +  eqv  —  qc  —  qr) k.  Due  to  the  microphysical  processes  that  involve  phase 
change  in  the  water  content,  the  source/sink  term  S  at  the  right-hand  side  of  the  equations  of 
potential  temperature  and  water  tracers  must  be  computed.  These  processes  are  modeled  by  the 
Kessler  parameterization  [12].  Equations  (1)  must  be  solved  in  D  Vt  G  (0,  T).  Initial  and  boundary 
conditions  will  be  assigned.  8 ,  p,  and  p  are  related  through  the  equation  of  state  for  a  perfect  gas. 


2.1  Dynamic  dissipation  in  an  LES  sense 

The  boundedness  of  the  solution  computed  with  the  straightforward  SEM  approximation  of  (1) 
is  compromised  by  unphysical  Gibbs  oscillations.  To  stabilize  the  problem,  the  Euler  equations 
are  corrected  to  include  the  stresses  of  the  Navier-Stokes  equations,  where,  however,  the  viscosity 
coefficients  of  such  stresses  are  given  by  a  residual-based  approximation  that  leads  the  problem  to 
converge  to  the  entropy  solution,  as  proved  in  [16]. 

Remark  1  Because  a  saturation  adjustment  scheme  [19]  is  used  to  treat  the  moist  thermody¬ 
namics,  the  source  terms  are  set  to  zero  in  the  main  step  of  the  solution,  and  are  only  computed 
within  the  Kessler  sub-step.  For  this  reason,  the  sources  will  not  appear  in  the  regularized  version 
of  Equations  (1). 


We  write: 
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p't  +  u  •  Vp  +  pS7  •  u  = 
ut  +  u-Vu+^V-(I p)  = 

6t  +  u-ve  = 

q%t  +  u  ■Vqi  = 


V  •  {unV p) 

•  (pn  (Vu  +  VuT))  +  B 

V  •  (Kn\70) 

V  •  (i/cV0) . 


(2) 


Except  for  vc  that,  for  the  time  being,  is  set  to  a  constant,  the  viscosity  coefficients  that  appear  in 
the  first  five  equations  are  computed  dynamically  as  a  function  of  the  solution.  They  are  calculated 
element-wise  on  every  high  order  element  0e.  More  specifically,  given  the  sensible  temperature 
T  =  9(p/po)R/Cp  and  one  element  with  equivalent  length  JiQe .  we  start  by  defining  the  dynamic 
viscosities 


Mmax|ne  =  0.5/iQe  1 1  pn  I  loo,  He  |||un|  +  V7^"  ||oo,QeJ 

and 


(3) 


Ml|f2e 


hf)f  max 


P(p)llooA 

||pn-pn||oo,n 


P(u)lloo,ne  l|fi(g)lloo,ne  \ 

Un||oo,f2 '  ||0n-0n||oo,Qy 


(4) 


where  7  indicates  the  space  average  of  the  quantity  at  hand  over  and  the  ||  •  ||oo,fi  terms  at 
the  denominator  are  used  for  normalization  for  a  consistent  dimension  of  the  resulting  equation. 
Having  /imax  and  constructed,  we  can  compute  the  dynamics  coefficients  of  the  viscosity  terms 
in  Equations  (2)  as: 


Pnhe  =  min  (/j  max|f2e  i  Ml  |ne  )  >  ^n|ne  — 


Pr 

7-1 


Pr 


l|pn||oo,f2, 


P'nlcie  > 


(5) 


where  Pr  =  0.7  is  the  Prandtl  number  of  dry  air. 


Remark  2  To  keep  the  discussion  brief,  the  details  of  the  derivation  of  the  equations  is  not  re¬ 
ported  and  the  notation  is  somewhat  abused.  A  proper  formulation  will  be  reported  in  a  subsequent 
paper. 


3  Space  and  time  discretization 

Equations  (2)  are  approximated  in  space  by  high  order  spectral  elements  and  by  an  Implicit- 
Explicit  (IMEX)  method  in  time.  Details  can  be  found  in,  e.g.  [6]  (SEM)  and  [5]  (IMEX). 

4  Numerical  Tests 

The  SEM-LES  method  is  tested  against  benchmarks  of  ubiquitous  use  when  testing  new  atmo¬ 
spheric  dynamical  cores.  First,  the  model  is  verified  in  dry  mode.  We  perturb  a  neutrally  stable 
atmosphere  with  a  cold  thermal  anomaly  that  triggers  the  development  of  a  density  current.  Once 
we  have  verified  the  ability  of  the  model  to  handle  dry  dynamics,  we  solve  a  fully  three-dimensional 
supercell  triggered  by  the  thermal  perturbation  of  a  realistic,  moist,  partially  unstable  background 
state. 


4 


4000 


0  5000  10000  15000 

4000 


Figure  1:  Density  current:  9'  at  900  s.  Top:  128  x  1  x  32  el.  (A z  =  Ax  «  50  m).  Bottom:  256  x  1  x  64  el. 
(Az  =  Ax  «  25  m).  4t?l-order  elements. 


4.1  Density  current  in  a  pseudo-3D  domain 

The  density  current  is  a  standard  benchmark  in  the  development  of  atmospheric  codes  [20]. 
The  inviscid  version  of  [1]  is  used  for  our  analysis.  This  is  because  we  are  interested  in  assessing 
the  current  LES-like  approach  as  a  stabilizing  tool  that  does  not  require  further  viscosity.  The 
background  state  is  characterized  by  a  neutral  atmosphere  at  uniform  potential  temperature  0  =  300 
K  and  hydrostatically  balanced  pressure.  Due  to  the  symmetry  of  the  original  problem  with  respect 
to  the  plane  center  line  of  the  x  —  z  plane,  the  solution  is  computed  in  the  region  D  =  25.6  xoox 
6.4  km3  The  perturbation  6'  centered  in  (xc,xc)  =  (0,3)  km  has  radii  ( rx,ry,rz )  =  (4,  oo,2)km 
and  is  given  by  6'  =  0.5A#  (1  +  cos(ircR))  for  R  <  1,  with  amplitude  A 9  =  —15  K  and  section 
R  =  \/ (x  —  xc)/r%  +  (z  —  zc)/r%.  Periodic  boundary  conditions  are  used  along  y  whereas  no-flux 
conditions  are  set  in  x  and  z.  The  initial  velocity  is  zero  everywhere.  Figure  1  shows  the  fully 
developed  current  at  time  t  =  900  s  on  two  grids  with  uniform  resolutions  Ax  =  Az  =  50  m  and 
Ax  =  Ax  =  25  m.  To  measure  the  front  position  at  tf  =  900  s,  we  take  the  node  on  the  ground 
where  O'  =  —  1  K.  A  comparison  of  the  front  position  and  0'max  min  with  respect  to  previous  work  is 
reported  in  Table  1.  As  the  resolution  decreases,  the  front  appears  slower;  this  fact  is  also  observed 
in  Fig.  5  of  [20]. 

We  are  aiming  at  using  the  current  stabilizing  scheme  as  a  Large  Eddy  Simulation  scheme.  As 
a  first  analysis  in  this  direction,  we  compare  how  the  current  model  compares  with  the  classical 
model  by  Lilly  and  Smagorinsky  [13,  18].  The  Smagorinsky  solution  (implemented  within  NUMA 
as  well)  is  plotted  in  Figure  2.  A  more  thorough  and  quantitative  analysis  is  currently  being  carried 
out  by  the  authors.  At  a  resolution  Ax  =  Ax  ~  25m  and  by  plotting  comparable  contours  (values 
not  shown  in  the  plot),  the  two  models  are  highly  comparable,  although  the  degree  of  dissipation 
of  the  current  scheme  seems  lower  than  Smagorinsky’s  using  a  Smagorinsky  constant  Cs  =  0.14. 
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Figure  2:  Density  current  using  a  classical  Smagorinsky  SGS  scheme  with  constant  Cs  =  0.14:  6'  at  900s. 
256  x  1  x  64  el.  (Az  =  Ax  «  25m).  4th-order  elements. 


Significantly  more  sub-grid  structures  are  resolved  using  the  current  model.  Further  analysis  is 
though  required. 

Remark  3  Throughout  this  paper  we  have  discussed  an  LES  approach  to  stabilization.  Never¬ 
theless,  it  must  be  pointed  out  that  the  simulations  that  we  have  presented  are  not  necessarily  to 
be  viewed  as  LES  simulations  unless  finer  grids  are  used. 

4.2  3D  moist  convection 

The  three-dimensional  simulation  of  a  convective  cell  is  defined  in  the  domain  160  x  120  x  24  km3. 
The  initial  field  is  perturbed  by  a  temperature  anomaly  O'  3  K  warmer  than  the  surrounding 
environment,  which  is  given  by  the  sounding  of  [4],  The  domain  Qh  is  subdivided  into  40  x  30  x  24 
elements  of  order  4.  A  stretched  grid  along  z  is  used  to  make  the  resolution  higher  in  the  lower 
atmosphere  where  convection  is  triggered.  The  domain  is  crossed  by  a  horizontal  wind  along  the 
x-direction  with  a  12 ms-1  shear  at  z  =  2000m.  A  no-slip  condition  is  applied  on  the  surface 
boundary  while  periodic  boundaries  are  defined  along  x  and  y.  A  Rayleigh  type  absorbing  layer  is 
included  at  z  >  19000  m.  The  cloud  first  forms  at  approximately  500  s,  and  is  fully  develop  after 
4500  s.  A  3D  instantaneous  view  of  qc  is  plotted  in  Figure  3.  Qualitatively,  it  is  comparable  to 
previous  results  on  a  similar  case.  A  quantitative  evaluation  of  the  instantaneous  rain  on  the  ground 
is  plotted  in  Figure  4a,  whereas  the  cloud  content  obtained  by  averaging  qc  along  the  y— direction 
is  plotted  in  Figure  4b. 

5  Conclusions 

We  extended  to  high  order  spectral  elements  the  LES-based  stabilization  method  first  introduced 
in  [17]  for  the  finite  element  solution  of  fully  compressible  flows.  We  explored  the  capabilities  of 
this  inexpensive  technique  to  solve  the  Euler  equations  of  stratified  flows  at  the  low-Mach  regimes 
encountered  in  atmospheric  flows.  When  applied  to  dry  and  moist  simulations,  the  current  im¬ 
plementation  appears  to  give  satisfactory  results  that  are  comparable  to  others  presented  in  the 
literature.  Without  the  need  for  any  additional  viscosity,  this  dynamic  LES  scheme  proved  to  be 
sufficient  to  stabilize  the  spectral  element  solution  of  the  Euler  equations  in  atmospheric  applica¬ 
tions.  However,  since  a  thorough  analysis  was  not  carried  out  to  evaluate  this  approach  in  terms 
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tj  =  6000s.  qc. 


Figure  3:  Supercell:  3d  view  (using  az  =  -135  and  el  =  8)  of  qc  (grey  surface),  surface  velocity  (vectors), 
and  the  instantaneous  distribution  of  qr  on  the  ground  (contours). 
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(a)  Instantaneous  rain  distribution  on  the  ground  at 
t  =  6000  s. 


tf  =  6000s.  qCtave(g/kg). 
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(b)  Vertical  slice  of  the  distribution  of  qc  averaged 
along  the  y  direction. 


Figure  4:  3D  supercell:  horizontal  slice  of  qr  at  2  =  0  m  and  y-averaged  qc  at  t  =  6000  s. 
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Table  1:  Case  3.  Comparative  results  of  front  location  at  900s.  LES  (SEM),  VMS  (FE),  WRF-ARW 
V2.2  (FD),  f-wave  (FV),  filtered  Spectral  Elements  (SE),  filtered  Discontinuous  Galerkin  (DG),  REFC, 
REFQ  and  PPM  results  are  compared.  All  models  but  LES  and  VMS  used  artificial  diffusion  with  constant 
H  =  75  m2  s-1. 


Model 

Nel 

Order 

g  =  75  m2  s  1 

Front  Location  [m] 

LES  (25  m) 

256  x  1  x  64 

4*h 

NO 

15080 

LES  (50  m) 

128  x  1  x  32 

4*h 

NO 

14888 

LES  (100  m) 

64  x  1  x  16 

4th 

NO 

14546 

LES  (200  m) 

32  x  1  x  8 

4*h 

NO 

13736 

LES 

32  x  1  x  8 

6th 

NO 

14568 

LES 

32  x  1  x  8 

8th 

NO 

14754 

VMS  [15]  (25  m) 

NO 

14890 

VMS  [15]  (50  m) 

NO 

14629 

VMS  [15]  (75  m) 

NO 

14487 

VMS  [15]  (100  m) 

NO 

14355 

WRF-ARW  50  m 

YES 

14470 

SE  [6]  50m 

YES 

14767 

DG  [6]  50m 

YES 

14767 

f-wave  (FV)  [1]  50  m 

YES 

14975 

REFC  [20]  50  m 

YES 

14437 

PPM  [20]  50  m 

YES 

15027 

of  its  turbulence  modeling  properties,  much  additional  work  is  necessary  to  fully  assess  it  in  its 
applicability  as  a  turbulence  closure  for  atmospheric  simulations. 
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